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Critical states are sometimes identified experimentally through power-law statistics or universal 
scaling functions. We show here that such features naturally emerge from networks in self-sustained 
irregular regimes away from criticality. In these regimes, statistical physics theory of large interacting 
systems predict a regime where the nodes have independent and identically distributed dynamics. 
We thus investigated the statistics of a system in which units are replaced by independent stochastic 
surrogates, and found the same power-law statistics, indicating that these are not sufficient to 
establish criticality. We rather suggest that these are universal features of large-scale networks 
when considered macroscopically. These results put caution on the interpretation of scaling laws 
found in nature. 

PACS numbers: 02.50.-r, 05.40.-a, 87.18.Tt, 87.19.11, 87.19.1c, 87.19.1m 


I. INTRODUCTION 

Power law statistics are ubiquitous in physics and life 
sciences. They are observed in a wide range of com¬ 
plex phenomena arising in natural systems, from sandpile 
avalanches to forest fires mm, earthquakes amplitude so¬ 
lar flares, website frequent at ion as well as economics |3]. 
These distributions reveal unusual properties: the ob¬ 
served quantity has no typical scale, and is not well char¬ 
acterized by its mean (when it exists): it is distributed 
over orders of magnitudes and large deviations are not 
exceptionally rare, in the sense that extreme events are 
far more likely than they would be, for instance, in a 
Gaussian distribution. These singular properties, com¬ 
bined with their ubiquity in nature, has attracted wide 
attention in applied science. 

A number of theories were proposed in order to ac¬ 
count for the presence of such power-law distributions. 
Some theories use an analogy with statistical physics sys¬ 
tems and associate the presence of power-law scalings as 
the hallmark that the system could be operating at a 
phase transition. These theories associate power-law to 
a so-called notion of criticality : classically, in statisti¬ 
cal physics, critical phenomena are the behaviors occur¬ 
ring in systems in association with second order phase 
transitions. These are largely thought to be universal 
although no proof was provided yet. Indeed, a num¬ 
ber of statistics are found in vastly distinct models, and 
particularly in the Ising model [3] poised at its phase 
transition. The critical regime thus occurs only at very 
specific parameters values. At this regime the proper¬ 
ties of the system are particularly singular; in particular, 
a number of statistics are scale-invariant including for 
instance the size, duration and shape of collective phe¬ 
nomena. The seminal work of Bak, Tang and Wiesen- 
feld on the Abelian sandpile model [5] largely popular¬ 
ized the hypothesis that criticality may be the origin of 
power-laws observed in nature. Indeed, despite the fact 
that parameters for criticality are very rare, systems may 
self-organize naturally at criticality (at the phase transi¬ 
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tion point) without requiring fine tuning by a mechanism 
called self-organized criticality (SOC) [lj. This remark¬ 
able theory has sometimes led to the conclusion that nat¬ 
ural systems were critical based on the identification of 
power-law relationships in empirical data (see [2j|5] for 
reviews). 

While power laws, as well as phase transition (thus crit¬ 
icality) are well identified in models, a number of authors 
have underlined the importance of being cautious when 
claiming power-law behavior in finite systems, question¬ 
ing their relevance or usefulness mm- In particular, 
Stumpf and Porter in [bj aptly noted the importance to 
take a nuanced approach between theoretical and empir¬ 
ical statistical support reporting a power-law, as theories 
arise from infinite systems while and real systems and 
usual datasets are finite. 

In physics, several alternative theories were proposed 
to account for the presence of power-laws (see e.g. [3] for 
a review). In particular, it was noticed very early that 
a pure random typewriter (a monkey sitting at a type¬ 
writer) would generate texts with a power-law distribu¬ 
tion of word frequencies [8- identical as the one observed 
in data. This work brilliantly showed that power-law 
may arise from purely stochastic mechanisms, evidencing 
that some power-laws distributions may not reflect deep 
structures in the data. Li [9] formalized Miller’s theory, 
highlighting the fact that combinations of exponentials 
naturally yield power-law relationships. Bouchaud and 
others [TQlfTT] noted that inverse of regularly distributed 
quantities may show power-law distributions. Newman 
and colleagues showed that random walks generate sev¬ 
eral statistics scaling as power-laws [3], Yule introduced 
a process with broad applications, particularly in evolu¬ 
tion, naturally associated to power-law distributions, and 
Takayasu and collaborators mm showed that systems 
with aggregation and injection naturally generate clus¬ 
ters of size scaling as a power-law with exponent —3/2. In 
neuroscience, Benayoun, Wallace and Cowan have shown 
that neuronal networks models in a regime of balance of 
excitation and inhibition also provide power-law scalings 
of avalanches m- All these mechanisms are independent 
of any phase transition and arise away from criticality 
from a particular way of considering a random process. 

The hypothesis that networks of the brain operate at 
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criticality was introduced a decade ago with the develop¬ 
ment of recording techniques of local populations of cells 
and the analysis of specific events corresponding to col¬ 
lective bursts of activity separated by periods of silence. 
The first empirical evidence that neuronal avalanches 
may show power-law distributions of duration or size was 
derived from the analysis of neuronal cultures in vitro ac¬ 
tivity m- Based on an analogy with the sandpile model, 
these bursts were seen as “neuronal avalanches”, and 
were apparently distributed as a power-law with a slope 
close to —3/2, consistent with the distribution of sandpile 
avalanches. These in vitro findings were based on indi¬ 
rect evidences of spiking derived from local field poten¬ 
tials, extracellular signals associated with the summation 
of postsynaptic potentials (bursts produce negative peaks 
in the LFP signals) and affected by a number of events 
unrelated with spiking activity. Similar LFP statistics 
were later found in vivo in the awake monkey m and in 
the anesthetized rat |18j . These empirical evidences were 
used to draw strong conclusions on neural coding: the 
presence such power-laws would ensure maximized infor¬ 
mation capacity, dynamic range and information trans¬ 
mission HU [20]. However, the method of analyzing the 
amplitude of negative LFP peaks was shown to produce 
spurious power laws scalings m regardless of the spike 
activity of cells. Indeed, identical scalings were found in 
surrogate data, positive LFP peaks (that are indepen¬ 
dent of spiking activity), and also arise in elementary 
purely stochastic signals, such as excursions of Ornstein- 
Uhlenbeck processes through thresholds away from the 
mean, or in one-dimensional random walks [22] [23]: both 
duration and time of excursions show power-law statis¬ 
tics, and display shape invariance. It was further shown 
in that both in data and surrogate models, statistical 
significance of these power-laws of LFP peak was poor, 
and depended on the threshold chosen. In [24], Dehghani 
and collaborators have made a statistical analysis com¬ 
bining multielectrode in vivo recordings from the cerebral 
cortex of cat, monkey and human, and did not confirm 
the presence of power-laws. The data rather showed an 
optimal fit with two independent exponential processes. 

The poor statistical significance of LFP avalanche anal¬ 
ysis and the ambiguous results it yields has motivated an 
in-depth exploration of in vitro spiking data of cultures 
of neurons [25]. In this remarkable work, the authors 
used multi-unit data from in vitro cultures and addressed 
a number of properties of critical systems reported by 
Sethna et al unified theory of criticality of statistical 
systems [2U]. Of crucial importance in this theory are 
power-law scalings of specific events and the relationship 
between the different scaling exponents. Friedman and 
collaborators [25; revealed that the data was consistent 
with Sethna et al theory, since power-law scalings of both 
avalanche size and duration were reported, with slopes 
consistent with the critical exponents of —3/2 and —2 
respectively, but also the existence of a universal mean 
temporal profiles of firing events collapsing under specific 
scaling onto a single universal scaling function, thereby 
providing more substantial analogy between this in vitro 
system with statistical physics models at criticality. 

In the present paper, we show that these observations 
can arise naturally in neuronal systems that are not at 
criticality. We also provide a theoretical explanation for 
this, as well as analytic access to some of the relevant 
properties of such systems. 


II. AVALANCHES IN SPIKING NETWORK 
MODELS 

We start by investigating the avalanche distributions 
generated by the classicalmodel of spiking neuronal net¬ 
work with excitatory and inhibitory connections intro¬ 
duced by Nicolas Brunei in m- This model describes the 
interaction of n neurons described through their voltage 
(vi)i=i... n that decays to reversal potential in the absence 
of input, receives external input and spikes from intercon¬ 
nected cells, and fire action potentials when the voltage 
exceeds a threshold 0. In detail, the voltage of neuron i 
satisfies the equation: 

= ~ v i +- D ) (!) 
j =1 k >0 

while Vi < 0, and where r denotes the time constant 
of the membrane and R its resistivity. The input re¬ 
ceived by the neuron are Dirac (5s. Neuron i receives 
the kth impulse of neuron j, emitted at time £*■, after 
a delay denoted D and assumed constant, which alters 
its membrane potential of a quantity proportional to . 
Brunei’s model assumes that these coefficients are zero 
except for a fixed number of cells randomly chosen in 
the excitatory and inhibitory populations, for which the 
coefficient Jij have fixed values J and — gJ respectively 
(see m for details). This model is particularly inter¬ 
esting very for its versatility and ability and to produce 
diverse spiking patterns. One can classify the regimes of 
activity in terms of levels of synchrony and regularity, 
and different regimes emerge as a function of the relative 
levels of excitation and inhibition and the input, that can 
be identified through the computation of precise bifurca¬ 
tion curves m- The thus obtained regimes are termed 
activity states , to distinguish these from the statistical 
mechanics notion of phase : the regimes are here sepa¬ 
rated by bifurcations occurring in the mean-field limit of 
the system. Of special interest are the Asynchronous Ir¬ 
regular (Al) states, in which neurons fire in a Poisson-like 
manner, with no period of silence (hence no avalanche). 
This activity is evocative of the spike trains observed ex¬ 
perimentally in awake animals. Such sparsely-connected 
networks can also display periods of collective activity of 
broadly distributed duration interspersed by periods of 
silence, called Synchronous Irregular (SI) regimes, known 
to reproduce the qualitative features of spiking in anes¬ 
thetized animals or neuronal cultures. Although par¬ 
tially ordered and partially disordered, SI regimes occur 
for a wide range of parameter values (all in inhibition- 
dominated regimes) [27] and are not very sensitive to 
modifications of biophysical parameters. The SI regime 
is not at a transition in the activity regime; within this 
region, chaotic activity takes place. Indeed, inhibition 
dominates the excitation, thus when the activity spreads 
throughout the network, it triggers massive inhibition 
that naturally silences the network. This is what we ob¬ 
serve in simulations of the Brunei model (Fig. [l]) [28]. We 
investigate the statistics of spike units in both cases (see 
Fig. [if. 

First, in the Al state, the sustained and irregular firing 
does not leave room for repeated periods of quiescence at 
this network scale, preventing the definition of avalanches 
(see raster plot in Fig. 0B)) [2Q]. This sharply con¬ 
trasts with the SI state, in which we can define avalanches 
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FIG. 1. Avalanche spike statistics in the sparsely connected networks m with N = 5 000 neurons in the SI and AI states. 
(A) Raster plot of the SI together with the firing rate (below). (B) Raster plot in the AI state: spiking is asynchronous and 
faster (notice that the time window is shorter compared to A for legibility): no silence period arises. Parameters as in [23 
Fig.2D] with g — 4 (excitation-dominated regime) and v e xt/vthresh = 1.1. . In the SI states, (C) avalanche size (number 
of spikes) and (D) avalanche duration (number of bins) scale as power law (dashed line is the maximum likelihood fit), and 
averaged avalanche size scales as a power law with the duration avalanche duration according to the universal scaling law [26] 
(E). Average avalanche shapes collapse onto the same curve (F,G) very accurately. Avalanche shapes from 20 to 40 bins are 
plotted in gray, and 3 of the same size as in [25] are colored and with distinct thicknesses. Parameters as in [271 Fig. 2D] with 
N — 1000, g — 5 (inhibition-dominated regime) and Vext/vthresh — 0.9. 
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which display the same statistics assumed to reveal criti¬ 
cality in cultures [25]. Fig [IT A) represents the raster plot 
with typical avalanches taking place. Strikingly, both 
avalanche duration and avalanche size show excellent fit 
to a power law, validated by Kolmogorov Smirnov test of 
maximum likelihood estimator m, and the exponents 
found are consistent with those found in neuronal cul¬ 
tures [25]. In particular, we find that the size s of the 
avalanches (number of firings during one avalanche) ap¬ 
parently scales as s~ r with r = 1.42 (Fig. lJC)), close 
to the theoretical value of 3/2 (plotted on the figure for 
indication), and to the neuronal culture scaling (1.7 re¬ 
ported in [2S] ) • Using the Kolmogorov-Smirnov test m, 
we have tested the hypothesis that the data are dis¬ 
tributed as a power-law, and validated the hypothesis 
(Kolmogov-Smirnov distance 0.0097, p-value p > 0.99. 
This is also the case of the distribution of the duration 
t of the avalanches, found scaling as t~ a with a = 2.11 
(Fig.^D)), close to the theoretical value of 2 for critical 
systems and from the experimental value of 1.9 found in 
neuronal cultures [25]). The Kolmogorov-Smirnov dis¬ 
tance with a pure power law is very low, evaluated at 
0.031, which leads to a high p-value p = 0.99, validat¬ 
ing clearly the hypothesis of a power-law distribution of 
avalanche durations. Similar to what was observed ex¬ 
perimentally, the fit is valid on two octaves, and drops 
down beyond [31]. This exponential drop is classically re¬ 
lated to subsampling effects. We validated the presence 
of finite-size cutoffs by varying the size of the system, 
and found indeed that as the network size increases, the 
fit with a power-law distribution is valid on larger time 
intervals and the cutoff only arises at later times (see Sup¬ 
plementary Figure[ 6 |. Also consistent with neuronal data 
and crackling noise, we found that the average avalanche 
size scales very clearly as a power of its duration, but 
with a positive exponent 7 = 1.50, not consistent with 
the crackling noise relationship between exponents 


7 = 


a — 1 


r — 1 


( 2 ) 


that predicts an exponent equal to 7 = 2.64, but however 
quantitatively consistent with the in vitro data of m ■ 
We have also investigated the shape of the avalanches 
of different durations. We have found that, similarly to 
critical systems or to in vitro data, the avalanche shapes 
collapsed onto a universal scaling function (Fig. IT) when 
time is rescaled to a unit interval and shape rescaled by 
T 7_1 where 7 = 1.5 is the power law exponent of the 
average size. 

These scalings are not specific to the particular choice 
of parameters used in our simulations: we consistently 
find, in the whole range of parameter values correspond¬ 
ing to the SI state of Brunei’s model, and in particu¬ 
lar, away from any transitions, similar apparent power- 
law scalings with similar scaling exponents (see supple¬ 
mentary Figure [7|. Moreover, we confirmed, in addition 
to the fact that this regime is away from all transitions 
between the different network activity, that the system 
is not at criticality, by showing that relaxation towards 
the SI regime after perturbation is fast (within millisec¬ 
onds) within the region considered, although it does slow 
down close to the transition point (see supplementary 
Figure 8 |. We conclude that these statistics are valid 
in a whole regime of the system where the activity is 
synchronous irregular, and neither at a transition of the 


model nor in a regime consistent with the slow decay 
of perturbations associated to critical regimes. There¬ 
fore, finding power-law statistics in neuronal avalanches 
with exponents 1.5 and 2 do not reveal that the system 
operates at criticality, but rather seems a property of 
synchronous irregular states. 

The SI states are prominent in neuronal activity, espe¬ 
cially in anesthesia and neuronal cultures. It is precisely 
in these situations that power-law distributions of spike 
avalanches were reported experimentally [16] [18]. This 
regime differs from the awake activity where neurons fire 
in an AI manner. In these regimes, power-laws and crit¬ 
icality were reported based on LFP recordings [17]. We 
will come back to experimental evidences of power-laws 
in local-field potentials recordings of the activity in sec¬ 
tion (V] 


III. AVALANCHES AND BOLTZMANN 
MOLECULAR CHAOS PROPERTY 

The observation of such scaling relationship in simple 
models of neuronal networks away from criticality and 
for a broad range of parameter values suggests that the 
observed scaling is related to properties of the systems 
that are independent of the notion of criticality and that 
may be relatively general. These may therefore be related 
to the properties of the network activity, that we now 
describe in more detail. 


III.l. Propagation of chaos in neural networks 
models 

The classical theory of the thermodynamics of interact¬ 
ing particle systems states that in large networks (such 
as those of the brain), the correlations between neurons 
vanishes. This is also known as Boltzmann molecular 
chaos hypothesis ( Stofizahlansatz ) in reference to the hy¬ 
pothesis that the speed of distinct particles should be 
independent, key to the kinetic theory of gases of Lud¬ 
wig Boltzmann [32]. In mathematics, this property is 
called propagation of chaos , and is rigorously defined as 
follows: 

Definition 1 . For (Af , • • • , X™) a sequence of measures 
on (M d ) n . The sequence is said X-chaotic if for any k G 
N and ii, *-- , i*. a set of indexes independent of n, we 
have (X™ , • • • , X™ ) converge to k independent copies of 
X as n -A 00 . 

In our context, in the limit of large networks, neurons 
behave as independent jump processes with a common 
rate, which is solution of an implicit equation. This prop¬ 
erty is at the core of theoretical approaches to understand 
the dynamics of large-scale networks I22C2H35I. In the 
case of Brunei’s model, it is shown that since two neu¬ 
rons share a vanishing proportion of common input in 
the thermodynamic limit, allowing to conclude that the 
correlation of the fluctuating parts of the synaptic input 
of different neurons are negligible. This leads the authors 
to conclude that the spike trains of different neurons are 
independent point processes with an identical instanta¬ 
neous firing rate v(t) that is solution of a self-consistent 
equation (see [27] p. 186, first column]). In that view, ex¬ 
cept in the case of constant firing rate (the asynchronous 
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irregular state), neurons always show a certain degree of 
synchrony due to the correlations of the instantaneous 
firing rates of different neurons. 

Mathematically, several methods were developed for 
interacting particle systems and gases (see e.g. [36] for 
an excellent review). It is shown that generically, sys¬ 
tems of interacting agents, with sufficient regularity, show 
propagation of chaos. All these results are in particular 
valid for neuronal networks, as was shown recently in a 
number of distinct situations. Large n limits and prop¬ 
agation of chaos was demonstrated for large networks of 
integrate-and-fire neurons EH, firing-rate models with 
multiple populations [38], for conductance-based mod¬ 
els even in the large time regime 39], and was shown 
to hold in realistic network models incorporating delays 
and the spatial extension of the system [40J 4lj. Rigorous 
methods of convergence of particle systems show that the 
empirical measure of the system: 

1 n 

= n 52 Sxi 

3 = 1 

converges in law towards a unique solution. A very pow¬ 
erful and universal mathematical result demonstrated 
in m Lemma 3.1] ensures that that the convergence 
of the empirical measure of a particle systems towards 
a unique measure implies propagation of chaos. Several 
methods may be used to show convergence of the em¬ 
pirical measure, including in the case of neuroscience, 
coupling methods m , compacity estimates or large 
deviations [43] . 

These results indicate that a universal form of activity 
emerges from neural networks, whereby neurons are inde¬ 
pendent copies of the same process. Before investigating 
the avalanche statistics of such regimes of activity, let us 
discuss the plausibility of the existence of these regimes 
in neuronal data. 


III.2. Decorrelation in experimental data and 
models 

In natural environments, regularity of sensory input to 
the brain may create strong and long-range correlations 
in space and time mss. It soon appeared that these 
correlations would be detrimental for the brain to encode 
sensory stimuli and detect changes efficiently na mi. 
Theoretical models of the visual system in particular have 
shown that decrease of redundancy by decorrelation was 
important for efficient encoding of natural images E2B- 
[52] . This was confirmed experimentally. In [53], the 
authors used high density two-dimensional electrode ar¬ 
ray and found in particular a marked exponential decay 
of correlation of excitatory cells. An clear confirmation 
of decorrelation even for closeby cells receiving similar 
input was recently brought in a remarkable experimen¬ 
tal work where chronically implanted multielectrode ar¬ 
rays were developed and implanted in the visual cortex 
of the macaque [54]. This protocol produced exquisite 
data allowing to show that even nearby neurons, gen¬ 
erally thought to be strongly connected and to receive 
a substantial part of common input showed a very low 
level of correlation. Similar decorrelation results were re¬ 
ported in the rodent neocortex [5b with the same level 
of accuracy. 


The origin of this decorrelation is still controversial and 
several assumptions were formulated, including the role 
played by adaptation ionic currents that could play cen¬ 
tral role in temporal decorrelation [56], negative correla¬ 
tions associated with the co-evolution of excitatory and 
inhibitory cells activity m, or sophisticated and robust 
mechanisms relying on neuronal nonlinearities and ampli¬ 
fied by recurrent connectivities EH, that was compatible 
with pattern decorrelation observed in the olfactory bulb 
of the zebra fish. 

All these experimental findings confirm that regimes 
in which neurons are independent are plausible represen¬ 
tations of neural networks activity. We now investigate 
the avalanche statistics of such networks. 


III.3. Statistics of networks in the molecular chaos 
regime 

Both the mathematical analysis of neuronal network 
models and fine analysis of the structure of spike trains 
motivates the study of ensembles of neurons that are in¬ 
dependent but with common non-stationary statistics. 
The simplest model made of one could think of is to con¬ 
sider a collection of independent Poisson processes with 
identical time-dependent rates. 

In that view, cells with constant firing rates resemble 
AI regimes. To generate a stochastic surrogate of the 
SI regime, the common rate of the cells should display 
non-periodically periods of silence. An obvious choice 
would be to replay the rates extracted from the SI state, 
and indeed such a surrogate generated power-law statis¬ 
tics (not shown), but in this case we could not rule out 
whether the power-law statistics are encoded in the rate 
functions. To show that this is not the case, we gener¬ 
ated a surrogate independent of the rate functions, by 
using a common rate of firing of the neurons given by 
the positive part of the Ornstein-Uhlenbeck process: 

Pt = —&Pt + 

with (£ t ) a Gaussian white noise. This choice is inter¬ 
esting in that although periods of silence do not occur 
periodically, the duration between two such silences have 
a finite mean. Actually, the distribution of excursions 
shape and duration of the Ornstein-Uhlenbeck process 
are known in closed form [58]. These distributions are, 
of course, not heavy tailed: they have exponential tails 
with exponent a which is the timescale of decay of the 
process. 

We investigated the collective statistics of N = 2 000 
independent realizations of Poisson processes with this 
rate. The resulting raster plot is displayed in Fig. | 2 ][A). 
While the firing is an inhomogeneous Poisson process, 
macroscopic statistics showpower-law distributions for 
the size (r = 1.47, Fig. J2[B)) and for the duration 
(a = 1.9, Fig. [2|C)), both statistically significant [30] 
and consistent with critical exponents. Again, a linear 
linear relationship between average avalanche size and 
duration is found, with a coefficient evaluated to 7 = 1.4 
and explaining all the variance but 610 -4 . Again, this 
coefficient is not consistent with the crackling relation¬ 
ship ©• Notwithstanding, we found that the averaged 
shape of avalanches of a given duration collapse on one 
universal curve when the amplitude is rescaled by the 
duration to the power 7 — 1 (Fig. HD,E)). 
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time in avalanche t (ms) fraction of avalanche duration t/T 

FIG. 2. Avalanches statistics in the independent Poisson model with Ornstein-Uhlenbeck firing rate (a = a = 1). Apparent 
power-law scalings, together with scale invariance of avalanches shapes. Avalanches of size 10 to 40 in gray, 3 specific trajectories 
highlighted. 


Of course, these statistics are not related to the nature 
of the rate chosen. We display in Supplementary Fig. [9] 
for instance the avalanche statistics of independent Pois¬ 
son processes with rate given by the positive part of a 
reflected Brownian, and we find exactly the same power- 
law statistics of avalanche shapes and durations, as well 
as a very nice collapse of avalanche shapes. 


III.4. Analytical derivations in the stationary and 
slow rate regime 


where p denotes the joint probability for the firing rate 
to have a specific sequence of values (go, • • • , q T +i)- 
This formula remains quite complex. Assuming 
now that the rate is extremely slow compared to the 
avalanches, one can simplify further the probability of 
p® of an avalanche of size r: 

P°r = [ q 2 n (l~q n ) T Pl(q)dq, 

J 0 

thus with a simple change of variable: 


In this simple model, it is actually very simple to in¬ 
deed compute explicitly the distribution of avalanche du¬ 
ration for instance. Disregarding the Poisson nature of 
the firings and the type of rate chosen, we can indeed 
write down the probability for an avalanche of duration 
t to occur at a specific time. These distributions can be 
computed analytically in cases where the spiking is de¬ 
scribed by a point process. In that case, denoting by p(t) 
the probability for a neuron to spike in the time inter¬ 
val [£, £ + 5], the probability to observe an avalanche of 
size r starting at time t* in a collection of n independent 
realizations is 


1 

p® = — x 1+ ™(l - x) r pi(x™)dx. 

n Jo 

As expected, when n oc, the probability of having 
an avalanche of prescribed, finite duration goes to zero 
as 1/n. The typical shape of the distribution can be 
obtained by rescaling this probability by n. Since we are 
interested in the logarithmic shape of the distribution, 
we disregard any multiplicative constant. As n oo, 
we thus obtain that the probability profile converges 
towards a universal limit independent of the particular 
shape of the distribution p, precisely given by: 


q (t*r q (t*+T + l) n l[(l-q(t*+t) n ). 

t =1 



1 

(r + l)(r + 2) 


Assuming stationarity, the probability of finding an 
avalanche of duration r is given by 

[ QoQr+1 11(1 “ qt) d Pr+2{qO, ' , <7r+1) 

J[ o.i] T+2 t=l 


which is indeed a power law with exponent —2, identi¬ 
cal to the one arising in critical systems, consistent with 
those reported in neuronal data [25] . in neuron models 
(Fig. [TJ and in surrogate systems (Fig. [2|. 

We now show that this extends to the distribution of 
avalanche size and the scaling of the mean avalanche size 
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with duration. Indeed, the size of avalanches of duration 
r have a binomial distribution, corresponding to s — r 
successes among (n— l)r independent Bernoulli variables 
with probability of success l—q(t). Moreover, De Moivre- 
Laplace theorem m ensures convergence as n increases 
towards a normal variable with mean (n — 1 )rp(t) and 
variance (n — 1 )rp(t)q(t) with p(t) = 1 — q(t). Using 
again our separation of timescale and stationarity hy¬ 
potheses, we find the probability the probability of find¬ 
ing an avalanche of size s and duration r averaged on the 
firing rate: 


p n 


(s,r) - / 
Jo 


1 (s _ T _ n(1 _ <?)T) 2 q 2 nn_ q n\r 


which converges, as n —» oo, towards: 

^ (s-T(l + log(u))) 2 u(l — Uy 

^ Oh- lnD-f'n'l V J 


P°° 


(s,t) ~ f 
Jo 


g 2 t log(n) 


du 


v/ 27rr log(w) 

We thus obtain at leading order the size distribution: 


s gi 

) ~ e " e [ 

T—l J ° 


V(s 


{s-t) 2 u(Ju( 1 — u)) 7 


g 2T log(u) 


y/27TTl0g(u) 


dx. (3) 


It is hard to further simplify this formula, but it can 
be easily evaluated numerically. We depict the result of 
this computation in Fig. [3j and illustrate the apparent 
power-law scaling with slope —3/2. 



FIG. 3. Universal shape of the avalanche size distribution in 
the slow rate and large n limit. 


Eventually, we obtain for the average size A r of 
avalanches of duration r: 
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IV. SPIKE PATTERN ENTROPY AND 
INFORMATION CAPACITY 

We have thus proved that power-law distributions of 
avalanches do not necessarily reveal that the network is 
operating at criticality. However, a number of theories 
have proposed that operating at criticality was an opti¬ 
mal regime of information processing in the brain, maybe 
selected by evolution as a useful trait for the nervous sys¬ 
tem [I9l l60U64] . The question that arises is thus whether 
these theories break down when power-law statistics no 
more arise from the system operating at criticality, but 
from a mean-field Boltzmann chaos regime. 

In order to address this outstanding question, we came 
back to the methods used in order to demonstrate opti¬ 
mality of data processing capabilities at criticality. These 
theories rely on the computation of the information ca¬ 
pacity of the network in different activity regimes evalu¬ 
ated, following Shannon’s information theory, as the en¬ 
tropy of the patterns of spike fired. In detail, a spike 
pattern in a network of size N is a TV-uplet s G {0,1}^, 
with Si = 1 (resp. Si = 0) if neuron i has fired (resp., not 
fired) in a specific timebin. If p denotes the probability 
of occurrence of spike patterns, the entropy is given by: 

Entropy = E p(s) logO(s)). 

56(0,1}^ 

In order to test whether our theory accounting for the 
emergence of power-law distributions in the absence of 
criticality challenges high information capacity of neu¬ 
ronal networks, we computed the information capacity 
of Brunei’s model in different regimes (see Fig. [§. The 
numerical results show that this is not the case, and the 
information capacity is maximal in the SI regime where 
power-law statistics of avalanches were observed. How¬ 
ever, we observe no difference between entropy levels in 
the SI or AI states. Therefore, we conclude that the 
maximality of entropy is not necessarily related to the 
emergence of power-law statistics. 

These observations can be well understood heuristi- 
cally. Indeed, the entropy of spike patterns is a measure 
of the variability of possible spike patterns observed in 
the course of neuronal activity. The diagram of Fig. [4] is 
thus not surprising. Indeed, while the diversity of spike 
patterns is reduced in the highly synchronized regular 
regimes, it will be large in the irregular regimes, both 
synchronized asynchronous [65]. In other words, entropy 
is maximized within irregular regimes where more diverse 
patterns are fired, and this independently of the under¬ 
lying mechanisms supporting the emergence of the irreg¬ 
ular activity. 


V. LOCAL FIELD POTENTIALS ARE 
AMBIGUOUS MEASURES OF CRITICALITY 


We thus conclude that while a power-law relationship 
persists between A r and r, the scaling exponent is not 
related to the exponents of the power-law of size (3/2) 
and duration (2) distribution through Sethna’s crack¬ 
ling noise relationship, which would predict an exponent 
equal to 2. Importantly, we note that the exponent found 
here is quantitatively perfectly consistent with the expo¬ 
nent found in in vitro data [25] in the neural network 
model (Fig. [Tj) or surrogate Poisson system (Fig. |2|. 


We have thus shown that discarding the criticality as¬ 
sumption does not degrade the quantifications of network 
efficiency This being said, we are facing an apparent con¬ 
tradiction. Indeed, our theory provides an account for 
the presence of critical statistics in networks in the SI 
regime, but discards AI states as possibly having criti¬ 
cally distributed avalanches. Evidences of critical statis¬ 
tics of avalanches in-vivo in the awake brain have been 


















FIG. 4. Computation of the entropy (circles) of the network 
in the Brunei’s model m for distinct values of the input in¬ 
tensity v ex t and inhibition ratio g. Color code indicates the 
entropy amplitude. Parameters as in the original model [27] 
Fig. 2B]. The blue surface represents the boundary between 
SI and AI regimes, and the green surface separates SR regimes 
from SI regimes. We observe a clear transition from low en¬ 
tropy (SR) to high entropy (SI+AI). 


reported, but are more scarce and controversial. Unlike 
neuronal cultures, the activity in the awake brain does 
not display bursts separated by silences, but is sustained. 
Using a macroscopic measurement of neuronal activity, 
the Local Field Potential (LFP), power laws could be 
shown from the distribution of peaked events [TT] . The 
motivation to use negative LFP peaks to deduce infor¬ 
mation from the distribution of spike avalanches relied 
on the fact that the amplitude of these peaks correlated 
with firing activity mm- 

However, this monotonic relationship between the 
number of spikes and the number of spikes do not im¬ 
ply that there should exist a relationship between the 
distribution of peak amplitude and avalanches. More¬ 
over, it was shown that power-laws naturally emerge 
from the random nature of the signal and the threshold¬ 
ing procedure used in this analysis, and moreover these 
power-laws may not be statistically significant [21 . In¬ 
deed, no power-law scaling could be found from unit ac¬ 
tivity, which were better fit by double-exponential dis¬ 
tributions [24]. These analyses rather suggest that the 
power-law statistics of LFP peaks does not reflect scale- 
invariant neural activity. 

In an attempt to clarify this, we investigate here 
whether the distribution of LFP peaks can display power- 
law scaling in spiking networks or in their stochastic sur¬ 
rogates. To obtain a more biophysical model where LFP 
can be defined, we considered the current-based Vogels 
and Abbott model [66] which provides a biologically re¬ 
alistic model of spiking network displaying asynchronous 
irregular and synchronous regular states, and in which 
synaptic currents are described by exponentially decay¬ 
ing functions with excitation and inhibition having dis¬ 
tinct time constants (in place of Dirac impulses in the 
Brunei model 0 see [66] for details). Simulations of the 
model provide an instantaneous distribution of postsy- 
naptic currents, from which we computed LFP signals. 

In detail, we have considered a spatially extended neu¬ 
ral network of 5 000 units randomly located on a 2- 


dimensional square and satisfying the by the Vogels- 
Abbott model. We evaluated an LFP signal Vlfp 
from the postsynaptic currents according to Coulomb’s 

law [57] : 


Vlfp 


Re 
47r 



where Vlfp is the electric potential at the electrode po¬ 
sition, R e = 230 Dcm is the extracellular resistivity of 
brain tissue, Ij are the synaptic currents of and rj is the 
distance between Ij and electrode position. Remarkably, 
applying to this more sophisticate model the same pro¬ 
cedures as in the original paper m , we found that the 
method cannot distinguish between structured or non- 
structured activity: for a fixed firing rate and bin size, 
we have been comparing in Figure [5] the LFP statistics 
of the avalanche duration in a network of neurons in the 
AI regime (A) or in the SI regime (C) that shows partial 
order of the firing and saw no difference. We also com¬ 
pared the statistics to those of independent Poisson pro¬ 
cesses with constant rate (B): the three instances show 
power-law scaling of avalanche duration, with the same 
exponent, that seem to rather be related to the firing 
rate and bin size than to the form of the network activ¬ 
ity. The scaling coefficient is again, close from 3/2, and 
varies with bin size. Clearly, in this case, the dynam¬ 
ics of LFP peaks cannot distinguish between critical and 
non-critical regimes. 





Time t (ms) 



FIG. 5. Avalanche analysis defined from the macroscopic vari¬ 
able Vlfp of a network of integrate-and-fire neurons with 
exponential synapses. Networks displaying an SI (left) or AI 
(middle) state, and a purely stochastic surrogate (right), show 
similar macroscopic power-law LFP peak statistics with the 
same exponent close to -3/2. 


VI. DISCUSSION 

In this paper, we have evaluated how power-law statis¬ 
tics and universal scaling can arise in the absence of crit¬ 
icality. We first outline the novel contributions of the 
present manuscript, then we discuss their significance. 

The contributions of the present manuscript are the 
following: 

(i). We first investigated the avalanche statistics of 
spiking neural networks models (Section 0 - We 
showed that these networks display power-law 
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statistics with the critical exponents -3/2 and -2, as 
well as collapse of the shapes of avalanches. This 
observation is robust and valid for a wide region of 
parameters corresponding to synchronous irregular 
activity in such networks, and thus away from any 
transition point (i.e., away from criticality). 

(ii). This robust numerical observation led us to inves¬ 
tigate theoretically how such statistics can emerge 
from the collective activity of spiking networks. 
One ubiquitous property of large neural networks 
models is that neurons behave as independent pro¬ 
cesses with the same statistics. This property, 
termed Boltzmann’s molecular chaos regime in sta¬ 
tistical physics, or propagation of chaos in mathe¬ 
matics, is universal in the dynamics of large-scale 
networks. We reviewed such properties in models, 
as well as recent experimental evidence supporting 
this decorrelation between the dynamics of cells. 


(iii). We next tested the hypothesis that power-laws may 
emerge from such large systems of weakly corre¬ 
lated units with similar statistics. We introduced 
and investigated numerically the dynamics of a sur¬ 
rogate network made of indepe ndent neurons shar¬ 
ing the same statistics (Section|III.3). Surprisingly, 
despite the simplicity of these systems, they display 
precisely the same power-law statistics of avalanche 
duration and size with the same exponents and 
shape collapse as fully-connected networks. 


(iv). This surrogate model is simple enough to derive in 
closed-form the statistics of the duration and size of 
the avalanches. Using this analytic expression, we 
demonstrate that the ‘critical’ exponents emerge 
naturally in large-scale systems (operating within 
the Boltzmann chaos regi me), w ithout the need to 
invoke criticality (Section III.4| ). 


(v). Current literature indicate that criticality is an op¬ 
timal information processing regime for the brain. 
We considered the same measures of information 
capability (entropy of spike trains) on our biolog¬ 
ically plausible models and showed that indeed, 
information is maximized in the SI regime where 
power-laws emerge. However, we show that this is 
not exclusive of the SI regime but same levels of en¬ 
tropy are found in AI regimes where no avalanche 
can be defined. 


(vi). Finally, we also addressed the presence of power- 
law scalings observed in LFP recordings by sim¬ 
ulating such signals emerging from more realis¬ 
tic neuronal network models with excitatory and 
inhibitory cells. Surprisingly, we observed that 
power-laws with the same exponents as observed in 
experiments are found in all regimes tested. This 
exponents persisted when neurons were replaced by 
independent Poisson processes. This clearly indi¬ 
cates that power-law scalings in LFPs do not con¬ 
stitute any proof of criticality in the underlying sys¬ 
tem. 


All together, these numerical and theoretical findings 
provide a new interpretation for the emergence of power- 
law statistics in large-scale systems, independent of the 
notion of criticality. We propose to explain the emergence 
of such scaling based on Boltzmann molecular chaos 


regime, known to govern the dynamics of most large- 
scale interacting systems E3 ES \M- In other words, 
power-law and universal scaling functions can be due to 
a mean-field effect in systems made of a large number 
of interacting units. Of course, this theory does not hy¬ 
pothesize that the elements considered (here, neurons) 
are disconnected in reality. To the contrary, the fact 
that the critical exponents still resist the removal of in¬ 
terconnections shows that such exponents do not need 
criticality to be explained. 

The main mechanism explored here, Boltzmann’s 
molecular chaos, is a universal feature of many statistical 
systems. The very particular structure of different parti¬ 
cles activity it induces, namely statistical independence 
of the particles behavior together with a correlation in 
the law, may induce as we have observed, the same type 
of power-laws as in critical systems, with universal co¬ 
efficients that are consistent with those found in critical 
systems. In agreement with this theory, we have seen 
that the same statistics are reproduced by a sparsely con¬ 
nected network and a surrogate stochastic process where 
the periods of firing and silences are themselves gener¬ 
ated by another stochastic process. This interpretation 
suggests that similar scaling relationships shall arise in 
more realistic neural network models with fixed connec¬ 
tivity patterns, in particular including axonal propaga¬ 
tion delays (constant delays are already present in Brunei 
model), dendritic structure, spike frequency adaptation, 
non-instantaneous synaptic transmission. Indeed, most 
of these elements will make the intrinsic dynamics of each 
cell more complex, but we do not expect that this com¬ 
plexity could affect the fact that these systems operate 
within the Boltzmann molecular chaos regime. Notwith¬ 
standing, models including synaptic plasticity, which is 
the process by which the brain acquires skills and stores 
memories, may not belong to the class of systems de¬ 
scribed in this paper. Indeed, in such systems, the con¬ 
nectivity patterns vary depending on the pairwise corre¬ 
lations of cells activity, and this relationship may com¬ 
pete with the establishment of Boltzmann’s molecular 
chaos regime. While this may not occur in the adult 
brain where plasticity is much slower than neuronal ac¬ 
tivity, distinct phenomena not described by our model 
may occur occur during the critical periods of brain de¬ 
velopment, when plasticity occurs at a faster timescale. 
Further experimental and theoretical investigations are 
necessary to characterize avalanche distributions in these 
systems, as well as to compute correlation levels to test if 
the decorrelation characteristic of Boltzmann’s molecular 
chaos occurs. 

Interestingly, a network operating in Boltzmann’s 
molecular chaos regime can be interpreted as a high di¬ 
mensional system with hidden variables, as studied re¬ 
cently in [SBJESj. In these contributions, the authors in¬ 
vestigate the rank distribution of high dimensional data 
with hidden latent variable, and show that such sys¬ 
tems display Zipf law scaling (power-laws with slope —1 
in the rank distribution) that generically arise from en¬ 
tropy consideration, and using the elegant identity be¬ 
tween entropy and energy shown in PE]. While these de¬ 
velopments do not generalize here, Boltzmann molecular 
chaos provides a natural explanation for the emergence of 
weakly correlated units with similar probability laws: in 
the neural network system case, the common rate could 
be seen as a latent variable, and both independence and 
irregularity build up only from the interactions between 
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cells. As a result of this theory, revealing apparent power- 
law scaling with exponents of -3/2 and 2, as well as shape 
collapse, may be entirely explained statistically; in par¬ 
ticular, these criteria constitute no proof of criticality 
and experimental studies solely relying on them should 
be re-evaluated. 

This statement is even more true when it comes to 
macroscopic measurements such as the LFP: both SI 
and AI regimes, as well as stochastic surrogates, display 
power-law statistics in the distribution of LFP peaks. 
This shows that systems of weakly correlated units, 
or their stochastic surrogates, can generate power-law 
statistics when considered macroscopically. Here again, 
the power-law statistics tells nothing about the critical or 
non-critical nature of the underlying system. This poten¬ 
tially reconciles contradictory observations that macro¬ 
scopic brain variables display power-law scaling m 
while no sign of such power-law scaling was found in the 
units mm- More generally, these results also put cau¬ 
tion on the interpretation of power-law relations found in 
nature. 

A question that naturally emerges is how to distinguish 
power-laws due to criticality from those due to Boltz¬ 
mann’s molecular chaos regime. We used the previous 
observation [26] that a prominent characteristic of crit¬ 
icality beyond the presence of power-law scalings is the 
particular relationship one finds between the exponents. 
We found here that the power-law scalings emerging in 
the absence of criticality did not satisfy this relationship. 
We propose to use that criterion as a possible way to dis¬ 
tinguish between power-law scaling due to criticality or 
due to Boltzmann’s molecular chaos. 

In conclusion, we have shown here that stochastic mod¬ 
els can replicate many of the experimental observations 
about ’critical’ exponents, which demonstrates that not 
only power-law scaling is not enough to prove criticality, 
but that we need new and better methods to investi¬ 
gate this in experimental systems. The fact that that 
such exponents are seen for networks and for stochas¬ 
tic systems, shows that they apply to a large class of 
natural systems and may be more universal than previ¬ 
ously thought. As Georges Miller noticed in his seminal 
paper [ 8 ] , examining random text typed by virtual mon¬ 
keys, the texts produced may not be interesting, but have 
some of the statistical properties considered interesting 
when humans, rather than monkeys, hit the keys. Simi¬ 
larly, the present results show systems that can emulate 
the power-law scaling seen in brain activity, but with no 
criticality involved. We thus cannot conclude on whether 
the brain operates at criticality or not, but we need more 
elaborate methods to resolve this point. 
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Appendix A: Power law statistics and maximum 
likelihood fits for stationary data 

We review here the methods used to fit the power 
law distributions, that closely follow the methodology 
exposed in m This methodology applies to stationary 
data. Taking the logarithm of the probability density 
of a power-law random variable, we obtain log (p(x)) = 
—alog(x) + log(a). The histogram of the power-law 
therefore presents an affine relation in a log-log plot. 
For this reason, power-laws in empirical data are often 
studied by plotting the logarithm of the histogram as a 
function of the logarithm of the values of the random 
variable, and doing a linear regression to fit an affine line 
to through the data points (usually using a least-squares 
algorithm). This method dates back to Pareto in the 
19th century (see e.g. ESP- The evaluated point x m [ n 
corresponding to the point where the data start having a 
power-law distribution is mostly evaluated visually, but 
this method is very sensitive to noise (see e.g. [73] and 
references herein). The maximum likelihood estimator of 
the exponent parameter a corresponding to n data points 

dE ^min IS- 

ri _ 

a = \+n(Y j \ogf^ . 

v x min 7 

2=1 

The log-likelihood of the data for the estimated parame¬ 
ter value is: 

L{a\X) = n log -) - a log (■ 

The parameter x m i n is evaluated then by minimizing 
the Kolmogorov-Smirnov distance: 

KS = max \S{x) — P(x)\ 

where S(x) is the cumulative distribution function (CDF) 
of the data and P(x) is the CDF of the theoretical dis¬ 
tribution being fitted for the parameter that best fits the 
data for x > x m i n ), as proposed by Clauset and colleagues 
in [74] . In order to quantify the accuracy of the fit, we use 
a standard goodness-of-fit test which generates a p-value. 
This quantity characterizes the likelihood of obtaining a 
fit as good or better than that observed, if the hypothe¬ 
sized distribution is correct. This method involves sam¬ 
pling the fitted distribution to generate artificial data sets 
of size n, and then calculating the Kolmogorov-Smirnov 
distance between each data-set and the fitted distribu¬ 
tion, producing the distribution of Kolmogorov-Smirnov 
distances expected if the fitted distribution is the true 
distribution of the data. A p-value is then calculated as 
the proportion of artificial data showing a poorer fit than 
fitting the observed data set. When this value is close 
to 1 , the data set can be considered to be drawn from 
the fitted distribution, and if not, the hypothesis might 
be rejected. The smallest p-values often considered to 
validate the statistical test are taken between 0.1 and 
0.01. These values are computed following the method 
described in m, which in particular involves generating 
artificial samples through a Monte-Carlo procedure. 

These methods, very efficient for stationary data, fail 
to evaluate the tails of non-stationary data as is the case 
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of neuronal data. A weighted Kolmogorov-Smirnov test 
with a refined goodness of fit estimate valid up to extreme 
tails [76] . 


Appendix B: Subsampling effects 


Of course, any analysis of finite sequences of data is 
subject to subsampling effects. While these may be 
neglected for light-tailed data, they become prominent 
when it comes to assessing possible slow decay of the tails 
of a statistical sample distribution. These effects were 
discussed in detail in a number of contributions. In the 
context of neuronal avalanches, these effects were charac¬ 
terized in m and the results show indeed a modification 
of the slope with subsampling, together with exponential 
cutoffs pushed to larger sizes as sampling becomes finer. 





FIG. 7. Avalanche statistics for the Brunei model with ran¬ 
domly chosen parameters within the SI regime. 



FIG. 6. Subsampling effects. Statistics for a randomly ex¬ 
tracted subset of neurons of size n = 125, 250, 500 among 
1 000 neurons whose dynamics is described by Brunei’s model. 
Theoretical power-laws with critical exponents are displayed 
in black dashed lines. 


equation: 



-v' + RIiit) 


where A (t) is the total current reaching the soma at time 
t. These currents arrive from the synapses made with 
other cells within the network, as well as from connec¬ 
tions to neurons outside the network. It is assumed that 
each neuron receives C ext connections to from excitatory 
neurons outside the network, and that these synapses 
are activated by independent Poisson processes with rate 
v ext . The current receives by neuron i is thus the sum: 


We have confirmed these results in our own data. In 
Fig.© we have computed the distribution of avalanche 
size and duration when considering only a fraction of 
the neurons for the statistics. In detail, we have sim¬ 
ulated the Brunei model [27] with N = 1000. This 
yields a raster plot, from which we have extracted a ran¬ 
domly chosen subset of n neurons, with n = N/k for 
k G {2,4,8}. We indeed observed that an exponential 
cutoff is shifted towards larger sizes and slopes increase 
with the subsampling ratio k. 


Appendix C: Brunei’s model 

In our simulations, we have used the neuronal network 
model introduced by Brunei in m and have referred 
to the different dynamical regimes of this system. We 
review here the model, provide all parameters used in our 
simulations and show that the conclusions drawn in one 
example of the synchronous irregular (SI) state are valid 
for all parameters tested within this regime. The model 
describes the dynamics of N integrate-and-fire neurons, 
80% of which are excitatory and the others inhibitory. 
In the model, it is assumed that each neuron receives 
C = sN randomly chosen connections, that are assumed 
to uniformly arise from the excitatory and the inhibitory 
population, thus 80% of the incoming connections to any 
cell come from the excitatory population. The network 
is assumed to be sparsely connected, thus £ « 1. The 
depolarization v l of neuron i at the soma satisfies the 


C-\-C t 

Rh{t) = T J2 Jij 5(t — t k j — D) 

3=1 k 

where the sum is taken over all synapses, are the 
synaptic efficacies, tj are the spike times at synapse 
j of neuron i, and D is the typical transmission de¬ 
lay, considered homogeneous at all synapses for simplic¬ 
ity. In order to simplify further the model, it is as¬ 
sumed that = J > 0 for all excitatory synapses, and 
Jij = — g J < 0 for inhibitory synapses. The parameter g 
is relevant in that it controls the balance between excita¬ 
tion and inhibition: if g < 4, the network is dominated by 
excitation, and otherwise it is dominated by inhibition. 
The neuron i fires an action potential when v 1 reaches a 
fixed threshold #, and the depolarization of neuron i is 
instantaneously reset to a fixed value V r where it remains 
fixed during a refractory period r rp (during this period, 
the neuron is insensitive to any stimulation). An impor¬ 
tant parameter is the ratio between the rate of external 
input v ext and the quantity denoted Vthresh correspond¬ 
ing to the minimal frequency that can drive one neuron, 
disconnected from the network, to fire an action poten¬ 
tial: Vthresh = oTJt (th e coefficient 0.8 in that formula 
corresponds to the fraction of excitatory neurons). 

In this model, the parameters that are kept free are 
the balance between excitation and inhibition g and the 
external firing rate v ext . All other parameters are chosen 
as in table [I] Using a mean-field analysis together with a 
diffusion approximation, the authors find that all neurons 
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£ 

D 

J 

T r p 

0 

V r 

0.1 

1.8 ms 

0.2 mV 

2 ms 

20 mV 

10 mV 


TABLE I. Para meters used in all simulations of Brunei’s 
model, as in 1221 - 


are independent point process driven by a common rate 
v(t) given by a self-consistent equation. Heuristically, 
during the time interval [t,t + dt\, the probability for any 
given to spike is given by v(t) dt , and the realization of 
this random variable are independent in the different neu¬ 
rons. When the rate v(t) depends on time, neurons thus 
show a level of synchrony, and when v(t) is constant, the 
regime is called asynchronous, in the parlance of [27] , In 
that paper, an analysis of the self-consistent rate equa¬ 
tion in the mean field limit led to the identification of 
several regimes that are depicted in Fig. [7] 

• The asynchronous irregular (AI) state in which v(t) 
converges towards a strictly positive constant value, 
which occurs when the excitatory external inputs 
are sufficiently large (is ext > vthresh ) and when in¬ 
hibition dominates excitation. 

• The synchronous regular (SR) regime corresponds 
to a state in which v(t) is a periodic function of 
time. This regimes arises in the excitation domi¬ 
nated regime, and the oscillations frequency is con¬ 
trolled essentially by the transmission delay D and 
the refractory period r rp (approximately varying as 
r rp /D). The transition thus occur close from the 
line g = 4. 

• The synchronous irregular (SI) regime occurs es¬ 
sentially in the inhibition-dominated regime when 
the input are not sufficient to drive the network to 
a sustained firing state, i.e. when v ext < ^thresh- 

We have reproduced in Fig. [7]the bifurcation diagram m 
Fig.2 B] with the bifurcation lines between AI, SI and SR 
states. Within the SI state, we have been randomly draw¬ 
ing 30 parameter points and analyzed the avalanches aris¬ 
ing for these parameters. We have found that all regimes 
show a very clear power-law distribution of avalanche size 
and duration with exponents consistent with the expo¬ 
nents — 1.5 and —2 predicted by the theory. 


Appendix D: No slowing down within the SI regime 

In addition to the fact that the SI regime is away from 
any transition between the different network regimes, we 
confirmed that the system did not show the typical prop¬ 
erties of critical states. A number of criteria were pro¬ 
posed in the Ising model to be characteristic of the critical 
regime. These include the divergence of the correlation 
length, of the heat capacity or magnetic susceptibility, 
that are all related to long-range correlations between 
spins. Here, the absence of order parameter and spatial 
dimension prevents from using similar criteria to inves¬ 
tigate the presence of critical dynamics. However, a cri¬ 
terion independent of the definition of an analogous of 
order parameter or magnetic susceptibility is the critical 


slowing-down occurring at phase transitions for dynami¬ 
cal systems. This criterion states that the relaxation time 
of the system, namely the time it takes for the system to 
return to its stationary regime after a perturbation, di¬ 
verges at criticality. 

In the present case, computing relaxation times is a 
challenge since the system is not at an equilibrium but 
within a chaotic regime, thus all perturbations produce 
massive changes in the dynamics of the system. Follow¬ 
ing the methodology developed in mm, we designed a 
numerical criterion to evaluate relaxation time to the SI 
regime. That regime is essentially defined by the alterna¬ 
tion of periods of collective activity followed by silences. 
We have thus perturbed the system by adding a con¬ 
stant input within a short time window (See Fig. [8FA)), 
which has the effect of switching the system into an 
asynchronous irregular regime where the firing is unin¬ 
terrupted. As the perturbation stops, the system quickly 
returns to an SI regime with alternations of silences and 
collective bursts. An upper bound of the relaxation time 
can thus be defined as the first time, after the pertur¬ 
bation has stopped, at which the system is completely 
silent. We have made extensive simulations within the 
SI regime to compute the relaxation time and obtained 
that the system returns to SI statistics after a few mil¬ 
liseconds after stimulation (on the order of 2ms). This 
time increases very fast close to the SI-AI transition as 
expected from the theory, but within the SI regime, the 
system did not show any indication of critical slowing- 
down. 


Appendix E: Diverse regimes of independent 
processes 

We have confirmed that the statistics of independent 
Poisson processes with fluctuating instantaneous firing 
rates produce avalanches with power-law distributions of 
avalanche size and durations, consistent with our theory. 
To this purpose, we have performed a similar analysis as 
in Fig. [^replacing Ornstein-Uhlenbeck firing rates by the 
positive part of a Brownian motion reflected at ±1. This 
choice was motivated by two constraints: the positive 
part was taken in order to consider only positive firing 
rates for consistency, and the reflection at ±1 was forced 
in order to prevent from having to long excursions of the 
Brownian motion, so that we can indeed assess that the 
heavy tails of the avalanche distributions are rather due 
to the statistical structure of the firings rather than due 
to possible very long excursions of the Brownian motion. 
The results of the simulations are provided in Fig. M As 
in the case of the positive part of the Ornstein Uhlen- 
beck process, we find very clear power-law distributions 
of avalanche size and durations, with slopes consistent 
with our theory, and a very clear collapse of the avalanche 
shapes. 

We add that beyond the collapse of the avalanche 
trajectories, the shapes onto these avalanche collapse 
may convey important information, as noted and in¬ 
vestigated in the context of one-dimensional random 
walks [22] [23]. We observe indeed that the shape of 
the network-generated avalanches are not similar to the 
shapes obtained in the Brownian or Ornstein-Uhlenbeck 
case, and may similarly contain an information that goes 
beyond pure shape collapse reported in neural data [25|. 
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FIG. 8. Relaxation time of the Brunei model in the SI regime. (A) Typical trajectory (blue) perturbed by a constant input 
(red) into a SI regime, returns to the SI state after a few milliseconds. This is true over the whole SI domain (B): relaxation 
times are on the order of a few milliseconds, and increase sharply close from the transition. 







FIG. 9. Avalanche statistics and shape collapses for indepen¬ 
dent Poisson processes with rates given by a reflected Brow¬ 
nian motion. 
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